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Abstract. We present a Monte Carlo study of the helicity moduli of an anisotropic 
classical three-dimensional (3D) XY-model of YBCO in superconducting state. It is 
found that both the a6-plane and the c-axis helicity moduli, which are proportional 
to the inverse square of the corresponding magnetic field penetration depth, vary 
linearly with temperature at low temperatures. The result for the c-axis helicity 
modulus is in disagreement with the experiments on high quality samples of YBCO. 
Thus we conclude that purely classical phase fluctuations of the superconducting order 
parameter cannot account for the observed c-axis electrodynamics of YBCO. 
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1. Introduction 

Emery and Kivelson presented [1,2] very strong physical arguments that in systems with 
low superfTuid density n s (e.g., high-T c copper oxide superconductors and the organic 
superconductors) the fluctuations in phase of the superconducting order parameter play 
a significant role. In essence, the superfluid density is proportional to the helicity mod- 
ulus of a superconductor [3], which measures the stiffness of the superconductor with 
respect to twists in phase of the order parameter; hence a low value of n s implies a 
low value of the maximum possible phase ordering temperature T^. If is compa- 
rable to the measured physical superconducting temperature T c , the observed super- 
conducting properties may be very different from what is predicted by the mean-field 
BCS/Eliashberg theory. They further argued that in the case of poor screening, as indi- 
cated by a low conductivity, the Coulomb interaction suppresses the local Cooper pair 
density fluctuations An s , which in turn enhances the phase fluctuations A(f> (An s A0 > 
1/2). Thus, the phase fluctuations might influence the superconducting properties over 
a wide range of temperatures below and above T c . Indeed, measurements of the a6-plane 
magnetic field penetration depth A(T) (A~ 2 oc n s ) on YBa2Cu306.95 [4] gave unambigu- 
ous evidence for three-dimensional (3D) XY critical scaling behavior in a temperature 
interval of about 10 K below T c =92.74 K. Roddick and Stroud [5] were the first to 
point out that classical phase fluctuations in a nodeless order parameter could produce 
a low-temperature A(T) which varies linearly with T in an isotropic three-dimensional 
superconductor. They also showed that this dependence persists when combined effects 
associated with Coulomb energy of local charge density fluctuations and Ohmic dissipa- 
tion are included in the model. This important result showed that the observed linear 
temperature dependence of the a6-plane A~ 2 in YBCO at low T [6] does not necessarily 
originate from quasiparticle excitations near the lines of the gap nodes at the Fermi 
surface. The same idea was advanced independently by Emery and Kivelson [1], who 
also argued that a large a6-plane conductivity in YBCO implies that the phase fluctua- 
tions are predominantly classical down to low temperature. More recently, Emery and 
Kivelson [7] argued that very extensive angle-resolved photoemission spectroscopy ex- 
periments found no evidence of quasiparticle excitations near the nodes of the gap even 
just below T c , and therefore the existence of nodes is not responsible for the observed 
linear T-dependence of the afr-plane A -2 in YBCO at low temperatures. 

If, indeed, the phase fluctuations in the order parameter are responsible for the 
observed linear temperature dependence of the a6-plane superfluid density at low T 
[1,5,7], then the experiments [8] on c-axis electrodynamics in YBCO are puzzling and 
have to be explained. Namely, it was found that the c-axis penetration depth A C (T) never 
has the linear temperature dependence observed in afr-plane. The best fit to the data up 
to about 40 K gives AA C (T) oc T 2 - 1 . If the phase fluctuations are predominantly classical 
down to low temperature one would expect the superfluid density to be linear at low T 
in all directions. Indeed, that is exactly what was obtained by Roddick and Stroud [5] 
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in a simple model for an anisotropic superconductor using a variational self-consistent 
phase phonon (SCPP) approximation. Here we present a detailed Monte Carlo (MC) 
study of helicity moduli for an anisotropic 3D XY-model which incorporates the bilayer 
structure of YBa 2 Cu307. As expected, we find the helicity moduli to decrease linearly at 
low temperatures in all directions. Thus the classical anisotropic 3D XY-model cannot 
account for the observed temperature dependence of the magnetic field penetration 
depth of YBa 2 Cu 3 7 along the c-axis. 

The rest of the paper is organized as follows. In section 2 we discuss our model 
and some details of the MC calculation. Section 3 contains our numerical results for the 
helicity moduli and our conclusions. 

2. The model and the simulation 

We consider an anisotropic 3D XY-model for a system with bilayer structure described 
by the Hamiltonian 

H =H ab + H c , (1) 

H ab = c °s(0M,s - <t>j,i,s) - h cos(2(0 MjS - <f>j,l,s))] , (2) 

l,s,{h3) 

H c = l~ J -L COS (0M,2 - 4>i,l,l) - J '± COs(0^ +:Li i - 4>i,l,2)\ ■ ( 3 ) 
l,i 

Here, the sum over I runs over a stack of bilayers, the sum over s = 1,2 runs over 
two layers in a given bilayer, denotes the nearest neighbors within a single layer, 
the sum over % runs over the sites in a given layer and <pi t i tS is the phase of the order 
parameter on site % of the layer s in the bilayer /. The ab-plane Josephson couplings with 
constants J\ and J 2 correspond to the transfer of one and two Cooper pairs [9] , Jj_ is the 
Josephson coupling constant between two layers within a given bilayer, and J' ± is the 
Josephson coupling constant between layers in two adjacent bilayers. This Hamiltonian 
is obtained from the coarse-grained Ginzburg-Landau model for YBCO [10] by ignoring 
the temperature dependence of the amplitude of the order parameter. 

The computation of the helicity modulus along direction perpendicular to the 
bilayers requires care as the separation between the layers in a bilayer, q,, is, in general, 
different from the separation between the bilayers, c' . If a uniform vector potential 
A is applied its effect on the Hamiltonian is to shift the phase difference between 
points ri and r 2 by A 1>2 = 2nA ■ (r 2 — r 1 )/$ ) where <3> — hc/2e is the flux quantum: 
0ri — 0r 2 — > — 0r 2 + ^.i,2- The momentum conjugate to <f>ij tS is the charge n^ iS 
(in the units of 2e) and its rate of change, given by —dH/d^i^, is the total current 
flowing into the site (i, I, s). The resulting expression can be used to define the current 
along each bond connected to site (i,l,s). In equilibrium (hi t i tS ) = and one obtains 
the average current conservation at each site (i, I, s). For a uniform A perpendicular to 
the layers, the average currents parallel to the bilayers vanish by symmetry and one has 

27T 2ty 

(J_l sinO^ - + -T-Ac b )) = (J' ± sin^j+i,! - i>Jj2 + ^-Ac')) , (4) 
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where the angular brackets denote the statistical average for the Hamiltonian H(A): 
(• • •) =j!id<p iU ■ ■■eM-PH(A))/Z(A), (3 = l/(k B T), Z{A) = j n# iiJ>a exp(-(3H(A)). 
The expression on the left hand side of Eq. (4) gives the average current from the site 
% in s — 1 layer of the bilayer I into the site % of s = 2 layer of the same bilayer. The 
expression on the right hand side of (4) gives the average current from the latter site 
into the site % of s — 1 layer in the adjacent layer l+l. 

Since the average currents along bonds parallel to A are equal, one can calculate 
the z-axis helicity modulus (the z-axis is perpendicular to the layers, i.e. it is along the 
crystallographic c-axis) by computing the derivative of any one of those average currents 
with respect to A at A = 0. For computational purposes using the Monte Carlo method 
we evaluate the average helicity modulus over all bonds parallel to A 

lzz = V Q \jp ^ [ J±c ^ cos (^.'. 2 ~ + ^i c '( c °s(0i,/+i,i - <t>i,iM 

~1zbT ^ Sill ^ i ''' 2 ~ + J>± SiU (^l+^ ~ ^.'.2)] 

x [^-LC 6 sin(0i i i i2 - 4>i,i,i) + J±c' sin^j+i,! - 4>i,i,2)\ ) 

l,i I 

+ k~^TJp ^ [ J ^ sin (^ 2 ~ ^.'.i)) + ■ 7 i( sin (& 1 i+i 1 i - (t>i,iM 

x [ J ± c &( sin (</W,2 - <Pi,i,i)) + J' ± c'{sm((f) ijl+1A - <f)i,i, 2 ))} \ ■ (5) 

l,i ) 

Here, N 3 is the number of bonds parallel to A. We note that computing ^y zz as d 2 F/dA 2 z) 
where F is the Helmholtz free energy, which is valid (up to a multiplicative constant) 
for orthorhombic lattices without a basis, would give a wrong result in our case. The 
in-plane helicity modulus j xx along x- (i. e. a-) direction is derived in analogous way. 

The elementary excitations of the XY-model are spin- waves and vortices [3], but in 
the zero temperature limit only the spin-waves are statistically relevant. It is useful to 
have the analytic expressions for ^ zz and ^ xx obtained from the spin-wave expansion in 
the zero temperature limit as a check on the Monte Carlo results at low temperatures. 
We find (in the units of (27t/^q)2N, where N is the number of tetragonal unit cells, and 
taking J 2 = for the sake of simplicity) 

I'M) = l - i(J ± c b + J' ± c') - (J± - J'i_)(J±c b - JW)/(J± + J'±)] - Oi z k h T , (6) 
where 



k 



( J ± c b + J' ± c') (— L- + 



(J ± QA(k) + JicV(k)) 



1 1 



v wi(k) w 2 (k). 



(7) 

«(k) = ( J ± + J' ± cos(k z c))/J(J ± + J' ± ) 2 - 4J ± J' ± sm 2 {k z c/2) , (8) 
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«'(k) = ( Jl cos(A; z c) + J' ± )/\f{J± + J' ± ) 2 - 4J ± Ji sin 2 (fc z c/2) , 



(9) 



Wl>2 (k) = 2 Jx sin 2 ^ + sin 2 ^ + ( J ± + Ji)/2 




T \\f{J± + J'±) 2 - 4J±A sin 2 (k z c/2) . 



(10) 



In Eq. (7) the sum over k runs over (k z > 0)-half of the Brillouin zone — ir < 
k x a, k y a, k z c < n and c = c& + d. The analytic expression for ^ xx is much simpler 
due to homogeneity of in-plane couplings and bond lengths 



To study the helicity moduli of the model described by the Hamiltonian (1-3) we 
have used Monte Carlo simulations based on the Metropolis algorithm [12]. In our 
finite-lattice MC simulations we have used periodic boundary conditions on lattices 
with 10x10x10, 20x20x20 and 30x30x30 sites. For a given set of parameters Ji, J 2 , 
Jl, q>, J' ± , d, the simulation would start at a low temperature with all phases aligned. 
The first 250 000 MC steps per site (sps) were thrown away followed by seven MC links 
of 250 000 MC sps each. At each temperature the range over which the phase angle was 
varied [13] was adjusted to ensure the MC acceptance rate of about 50%. The errors 
were calculated by breaking up each link into 5 blocks of 50 000 MC sps, then calculating 
the average values for each of 35 blocks and finally taking the standard deviation a of 
these 35 average values as an estimate of the error. The final configuration of the phase 
angles at a given temperature was used as a starting configuration for the next higher 
temperature. The cumulant analysis [13] was performed using the results for the three 
lattice sizes to determine the value of the transition temperature T c . 

3. Numerical results and conclusions 

We have performed MC simulations for different choices of parameters and since they 
give qualitatively similar results we present here only the data for a single set of 
parameters J\ = 1, J 2 = 0, J± = 0.1, J' ± = 0.04, q, = 0.4, d =1. The values of q, 
and d were chosen so that their ratio reflects the structure of YBa2Cus07_5 [14]. The 
value of Josephson coupling between the layers of a bilayer is taken to be one order 
of magnitude smaller than the in-plane Josephson coupling and the value of Josephson 
coupling between bilayers J' ± was chosen such that Jlq, = J' ± d. Since the helicity 
modulus is proportional to inverse square of the magnetic field penetration depth A(T) 
[5], we present our results for the helicity moduli in Fig. 1 as A 2 (0)/A 2 (T) as a function 
of T/T c . From the cumulant analysis of the numerical results for the three lattice sizes in 
in Fig. 1 we obtained ksT c = 1.265 (in the units of J±). This value should be compared 
to 0.89 - 0.95 for a two-dimensional square lattice [15-17] and 2.2 for a three-dimensional 
cubic lattice [18-20] (which we have also reproduced with our code) in the units of the 
coupling constant. 




(11) 
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Figure 1. The temperature dependence of helicity moduli for three different lattice 
sizes. 

As anticipated, A 2 (0)/A 2 (T) is decaying linearly at low T both along the a6-plane 
and along the c-axis with the slopes of -0.2366 ± 0.0001 and -0.415 ± 0.001, respectively, 
as determined by the linear fits to the low-temperature data (T <0.01T C ). These values 
should be compared with -0.16667 ± 0.00005 which we found for our isotropic 30x30x30 
lattice with the periodic boundary conditions. We note that the spin-wave expansion 
formulae (5-11) give the slopes of 0.2367 along afe-plane, -0.404 along c-axis for our 
anisotropic 3D XY-model, and -1/(6 Ji) = -0. 16666/ Ji for the isotropic 3D XY-model, 
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in excellent agreement with our Monte Carlo result. Note that our results for the 
afr-plane helicity modulus are somewhat lower than the spin-wave result 0.25/ Ji for 
two-dimensional square lattice, which is easily understood from Eqs. (10-11) when J± 
and J' ± are one to two orders of magnitude smaller than J\. 

In terms of the superfluid density fraction p s = A 2 (0)/A 2 (T) as a function of 
t = T/T c the slopes are -0.30 along afr-plane, -0.52 along c-axis for our anisotropic 
3D XY-model, and -0.37 for the isotropic 3D XY-model (compare with the value -0.4 
reported in [7]). 

Our calculations clearly demonstrate that classical phase fluctuations cannot 
explain the observed nonlinear temperature dependence of the c-axis magnetic field 
penetration depth for high-T c samples of YBCO [8]. 
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